The 2019-2020 Coronavirus Pandemic Analysis

Contact: Smith Research

BACKGROUND & APPROACH

I wanted to track and trend the coronavirus outbreak on my own curiosity. There are some interesting questions that may fall out of this, as it is a very historic moment, including scientifically and analytically (we have a large amount of data being shared across the globe, analyzed in real-time). The world has come to a halt because of it.
This analysis attempts to answer the following questions (more to come):

  1. What does the trend of the pandemic look like to date?
  2. What are future case predictions based on historical model?
  3. What interesting quirks or patterns emerge?

ASSUMPTIONS & LIMITATIONS: * This data is limited by the source. I realized early on that depending on source there were conflicting # of cases. Originally I was using JHU data… but this was always ‘ahead’ of the Our World In Data. I noticed that JHU’s website was buggy- you clicked on the U.S. stats but it didn’t reflect the U.S.. So I changed data sources to be more consistent with what is presented in the media (and Our World In Data has more extensive plots I can compare my own to). An interesting aside might be why the discrepancy? Was I missing something?
* Defintiions are important as is the idea that multiple varibales accumulate in things like total cases (more testing for example).

SOURCE RAW DATA: * https://ourworldindata.org/coronavirus
* https://github.com/CSSEGISandData/COVID-19/
*

INPUT DATA LOCATION: github (https://github.com/sbs87/coronavirus/tree/master/data)

OUTPUT DATA LOCATIOn: github (https://github.com/sbs87/coronavirus/tree/master/results)

TIMESTAMP

Start: ##—— Sun May 24 10:45:56 2020 ——##

PRE-ANALYSIS

The following sections are outside the scope of the ‘analysis’ but are still needed to prepare everything

UPSTREAM PROCESSING/ANALYSIS

  1. Google Mobility Scraping, script available at get_google_mobility.py
# Mobility data has to be extracted from Google PDF reports using a web scraping script (python , written by Peter Simone, https://github.com/petersim1/MIT_COVID19)

# See get_google_mobility.py for local script 

python3 get_google_mobility.py
# writes csv file of mobility data as "mobility.csv"

SET UP ENVIORNMENT

Load libraries and set global variables

# timestamp start
timestamp()
## ##------ Sun May 24 10:45:56 2020 ------##

# clear previous enviornment
rm(list = ls())

##------------------------------------------
## LIBRARIES
##------------------------------------------
library(plyr)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.0     ✓ purrr   0.3.3
## ✓ tibble  3.0.0     ✓ dplyr   0.8.5
## ✓ tidyr   1.0.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::arrange()   masks plyr::arrange()
## x purrr::compact()   masks plyr::compact()
## x dplyr::count()     masks plyr::count()
## x dplyr::failwith()  masks plyr::failwith()
## x dplyr::filter()    masks stats::filter()
## x dplyr::id()        masks plyr::id()
## x dplyr::lag()       masks stats::lag()
## x dplyr::mutate()    masks plyr::mutate()
## x dplyr::rename()    masks plyr::rename()
## x dplyr::summarise() masks plyr::summarise()
## x dplyr::summarize() masks plyr::summarize()
library(ggplot2)
library(reshape2)
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths
library(plot.utils)
library(utils)
library(knitr)

##------------------------------------------

##------------------------------------------
# GLOBAL VARIABLES
##------------------------------------------
user_name <- Sys.info()["user"]
working_dir <- paste0("/Users/", user_name, "/Projects/coronavirus/")  # don't forget trailing /
results_dir <- paste0(working_dir, "results/")  # assumes diretory exists
results_dir_custom <- paste0(results_dir, "custom/")  # assumes diretory exists


Corona_Cases.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_global.csv"
Corona_Cases.US.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_confirmed_US.csv"
Corona_Deaths.US.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_US.csv"
Corona_Deaths.source_url <- "https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_covid19_deaths_global.csv"

Corona_Cases.fn <- paste0(working_dir, "data/", basename(Corona_Cases.source_url))
Corona_Cases.US.fn <- paste0(working_dir, "data/", basename(Corona_Cases.US.source_url))
Corona_Deaths.fn <- paste0(working_dir, "data/", basename(Corona_Deaths.source_url))
Corona_Deaths.US.fn <- paste0(working_dir, "data/", basename(Corona_Deaths.US.source_url))
default_theme <- theme_bw() + theme(text = element_text(size = 14))  # fix this
##------------------------------------------

FUNCTIONS

List of functions

function_name description
prediction_model outputs case estumate for given log-linear moder parameters slope and intercept
make_long converts input data to long format (specialized cases)
name_overlaps outputs the column names intersection and set diffs of two data frame
find_linear_index finds the first date at which linearaity occurs
##------------------------------------------
## FUNCTION: prediction_model
##------------------------------------------
## --- //// ----
# Takes days vs log10 (case) linear model parameters and a set of days since 100 cases and outputs a dataframe with total number of predicted cases for those days
## --- //// ----
prediction_model<-function(m=1,b=0,days=1){
  total_cases<-m*days+b
  total_cases.log<-log(total_cases,10)
  prediction<-data.frame(days=days,Total_confirmed_cases_perstate=total_cases)
  return(prediction)
}
##------------------------------------------

##------------------------------------------
## FUNCTION: make_long
##------------------------------------------
## --- //// ----
# Takes wide-format case data and converts into long format, using date and total cases as variable/values. Also enforces standardization/assumes data struture naming by using fixed variable name, value name, id.vars, 
## --- //// ----
make_long<-function(data_in,variable.name = "Date",
                   value.name = "Total_confirmed_cases",
                   id.vars=c("case_type","Province.State","Country.Region","Lat","Long","City","Population")){

long_data<-melt(data_in,
                id.vars = id.vars,
                variable.name=variable.name,
                value.name=value.name)
return(long_data)

}
##------------------------------------------

## THIS WILL BE IN UTILS AT SOME POINT
name_overlaps<-function(df1,df2){
i<-intersect(names(df1),
names(df2))
sd1<-setdiff(names(df1),
names(df2))
sd2<-setdiff(names(df2),names(df1))
cat("intersection:\n",paste(i,"\n"))
cat("in df1 but not df2:\n",paste(sd1,"\n"))
cat("in df2 but not df1:\n",paste(sd2,"\n"))
return(list("int"=i,"sd_1_2"=sd1,"sd_2_1"=sd2))
}

##------------------------------------------

##------------------------------------------
## FUNCTION: find_linear_index
##------------------------------------------
## --- //// ----
# Find date at which total case data is linear (for a given data frame) 
## --- //// ----

find_linear_index<-function(tmp,running_avg=5){
  tmp$Total_confirmed_cases_perstate.log<-log(tmp$Total_confirmed_cases_perstate,2)
  derivative<-data.frame(matrix(nrow = nrow(tmp),ncol = 4))
  names(derivative)<-c("m.time","mm.time","cumsum","date")
  
  # First derivative
  for(t in 2:nrow(tmp)){
    slope.t<- tmp[t,"Total_confirmed_cases_perstate.log"]- tmp[t-1,"Total_confirmed_cases_perstate.log"]
    derivative[t,"m.time"]<-slope.t
    derivative[t,"date"]<-as.Date(tmp[t,"Date"])
  }
  
  # Second derivative
  for(t in 2:nrow(derivative)){
    slope.t<- derivative[t,"m.time"]- derivative[t-1,"m.time"]
    derivative[t,"mm.time"]<-slope.t
  }
  
  #Compute running sum of second derivative (window = 5). Choose point at which within 0.2
  for(t in running_avg:nrow(derivative)){
    slope.t<- sum(abs(derivative[t:(t-4),"mm.time"])<0.2,na.rm = T)
    derivative[t,"cumsum"]<-slope.t
  }
  
  #Find date -5 from the stablility point
  linear_begin<-min(derivative[!is.na(derivative$cumsum) & derivative$cumsum==running_avg,"date"])-running_avg
  
  return(linear_begin)
}

READ IN DATA

# Q: do we want to archive previous versions? Maybe an auto git mv?

##------------------------------------------
## Download and read in latest data from github
##------------------------------------------
download.file(Corona_Cases.source_url, destfile = Corona_Cases.fn)
Corona_Totals.raw <- read.csv(Corona_Cases.fn, header = T, stringsAsFactors = F)

download.file(Corona_Cases.US.source_url, destfile = Corona_Cases.US.fn)
Corona_Totals.US.raw <- read.csv(Corona_Cases.US.fn, header = T, stringsAsFactors = F)

download.file(Corona_Deaths.source_url, destfile = Corona_Deaths.fn)
Corona_Deaths.raw <- read.csv(Corona_Deaths.fn, header = T, stringsAsFactors = F)

download.file(Corona_Deaths.US.source_url, destfile = Corona_Deaths.US.fn)
Corona_Deaths.US.raw <- read.csv(Corona_Deaths.US.fn, header = T, stringsAsFactors = F)

# latest date on all data:
paste("US deaths:", names(Corona_Deaths.US.raw)[ncol(Corona_Deaths.US.raw)])
## [1] "US deaths: X5.23.20"
paste("US total:", names(Corona_Totals.US.raw)[ncol(Corona_Totals.US.raw)])
## [1] "US total: X5.23.20"
paste("World deaths:", names(Corona_Deaths.raw)[ncol(Corona_Deaths.raw)])
## [1] "World deaths: X5.23.20"
paste("World total:", names(Corona_Totals.raw)[ncol(Corona_Totals.raw)])
## [1] "World total: X5.23.20"

PROCESS DATA

  • Convert to long format
  • Fix date formatting/convert to numeric date
  • Log10 transform total # cases
##------------------------------------------
## Combine death and total data frames
##------------------------------------------
Corona_Totals.raw$case_type<-"total"
Corona_Totals.US.raw$case_type<-"total"
Corona_Deaths.raw$case_type<-"death"
Corona_Deaths.US.raw$case_type<-"death"

# for some reason, Population listed in US death file but not for other data... Weird. When combining, all datasets will have this column, but US deaths is the only useful one.  
Corona_Totals.US.raw$Population<-"NA" 
Corona_Totals.raw$Population<-"NA"
Corona_Deaths.raw$Population<-"NA"

Corona_Cases.raw<-rbind(Corona_Totals.raw,Corona_Deaths.raw)
Corona_Cases.US.raw<-rbind(Corona_Totals.US.raw,Corona_Deaths.US.raw)
#TODO: custom utils- setdiff, intersect names... option to output in merging too
##------------------------------------------
# prepare raw datasets for eventual combining
##------------------------------------------
Corona_Cases.raw$City<-"NA" # US-level data has Cities
Corona_Cases.US.raw$Country_Region<-"US_state" # To differentiate from World-level stats

Corona_Cases.US.raw<-plyr::rename(Corona_Cases.US.raw,c("Province_State"="Province.State",
                                                  "Country_Region"="Country.Region",
                                                  "Long_"="Long",
                                                  "Admin2"="City"))


##------------------------------------------
## Convert to long format
##------------------------------------------
#JHU has a gross file format. It's in wide format with each column is the date in MM/DD/YY. So read this in as raw data but trasnform it to be better suited for analysis
# Furthermore, the World and US level data is formatted differently, containing different columns, etc. Recitfy this and combine the world-level stats with U.S. level stats.

Corona_Cases.long<-rbind(make_long(select(Corona_Cases.US.raw,-c(UID,iso2,iso3,code3,FIPS,Combined_Key))),
make_long(Corona_Cases.raw))


##------------------------------------------
## Fix date formatting, convert to numeric date
##------------------------------------------
Corona_Cases.long$Date<-gsub(Corona_Cases.long$Date,pattern = "^X",replacement = "0") # leading 0 read in as X
Corona_Cases.long$Date<-gsub(Corona_Cases.long$Date,pattern = "20$",replacement = "2020") # ends in .20 and not 2020
Corona_Cases.long$Date<-as.Date(Corona_Cases.long$Date,format = "%m.%d.%y")
Corona_Cases.long$Date.numeric<-as.numeric(Corona_Cases.long$Date)

kable(table(select(Corona_Cases.long,c("Country.Region","case_type"))),caption = "Number of death and total case longitudinal datapoints per geographical region")
Number of death and total case longitudinal datapoints per geographical region
death total
Afghanistan 123 123
Albania 123 123
Algeria 123 123
Andorra 123 123
Angola 123 123
Antigua and Barbuda 123 123
Argentina 123 123
Armenia 123 123
Australia 984 984
Austria 123 123
Azerbaijan 123 123
Bahamas 123 123
Bahrain 123 123
Bangladesh 123 123
Barbados 123 123
Belarus 123 123
Belgium 123 123
Belize 123 123
Benin 123 123
Bhutan 123 123
Bolivia 123 123
Bosnia and Herzegovina 123 123
Botswana 123 123
Brazil 123 123
Brunei 123 123
Bulgaria 123 123
Burkina Faso 123 123
Burma 123 123
Burundi 123 123
Cabo Verde 123 123
Cambodia 123 123
Cameroon 123 123
Canada 1722 1722
Central African Republic 123 123
Chad 123 123
Chile 123 123
China 4059 4059
Colombia 123 123
Comoros 123 123
Congo (Brazzaville) 123 123
Congo (Kinshasa) 123 123
Costa Rica 123 123
Cote d’Ivoire 123 123
Croatia 123 123
Cuba 123 123
Cyprus 123 123
Czechia 123 123
Denmark 369 369
Diamond Princess 123 123
Djibouti 123 123
Dominica 123 123
Dominican Republic 123 123
Ecuador 123 123
Egypt 123 123
El Salvador 123 123
Equatorial Guinea 123 123
Eritrea 123 123
Estonia 123 123
Eswatini 123 123
Ethiopia 123 123
Fiji 123 123
Finland 123 123
France 1353 1353
Gabon 123 123
Gambia 123 123
Georgia 123 123
Germany 123 123
Ghana 123 123
Greece 123 123
Grenada 123 123
Guatemala 123 123
Guinea 123 123
Guinea-Bissau 123 123
Guyana 123 123
Haiti 123 123
Holy See 123 123
Honduras 123 123
Hungary 123 123
Iceland 123 123
India 123 123
Indonesia 123 123
Iran 123 123
Iraq 123 123
Ireland 123 123
Israel 123 123
Italy 123 123
Jamaica 123 123
Japan 123 123
Jordan 123 123
Kazakhstan 123 123
Kenya 123 123
Korea, South 123 123
Kosovo 123 123
Kuwait 123 123
Kyrgyzstan 123 123
Laos 123 123
Latvia 123 123
Lebanon 123 123
Lesotho 123 123
Liberia 123 123
Libya 123 123
Liechtenstein 123 123
Lithuania 123 123
Luxembourg 123 123
Madagascar 123 123
Malawi 123 123
Malaysia 123 123
Maldives 123 123
Mali 123 123
Malta 123 123
Mauritania 123 123
Mauritius 123 123
Mexico 123 123
Moldova 123 123
Monaco 123 123
Mongolia 123 123
Montenegro 123 123
Morocco 123 123
Mozambique 123 123
MS Zaandam 123 123
Namibia 123 123
Nepal 123 123
Netherlands 615 615
New Zealand 123 123
Nicaragua 123 123
Niger 123 123
Nigeria 123 123
North Macedonia 123 123
Norway 123 123
Oman 123 123
Pakistan 123 123
Panama 123 123
Papua New Guinea 123 123
Paraguay 123 123
Peru 123 123
Philippines 123 123
Poland 123 123
Portugal 123 123
Qatar 123 123
Romania 123 123
Russia 123 123
Rwanda 123 123
Saint Kitts and Nevis 123 123
Saint Lucia 123 123
Saint Vincent and the Grenadines 123 123
San Marino 123 123
Sao Tome and Principe 123 123
Saudi Arabia 123 123
Senegal 123 123
Serbia 123 123
Seychelles 123 123
Sierra Leone 123 123
Singapore 123 123
Slovakia 123 123
Slovenia 123 123
Somalia 123 123
South Africa 123 123
South Sudan 123 123
Spain 123 123
Sri Lanka 123 123
Sudan 123 123
Suriname 123 123
Sweden 123 123
Switzerland 123 123
Syria 123 123
Taiwan* 123 123
Tajikistan 123 123
Tanzania 123 123
Thailand 123 123
Timor-Leste 123 123
Togo 123 123
Trinidad and Tobago 123 123
Tunisia 123 123
Turkey 123 123
Uganda 123 123
Ukraine 123 123
United Arab Emirates 123 123
United Kingdom 1353 1353
Uruguay 123 123
US 123 123
US_state 401103 401103
Uzbekistan 123 123
Venezuela 123 123
Vietnam 123 123
West Bank and Gaza 123 123
Western Sahara 123 123
Yemen 123 123
Zambia 123 123
Zimbabwe 123 123
# Decouple population and lat/long data, refactor to make it more tidy
metadata_columns<-c("Lat","Long","Population")
metadata<-unique(select(filter(Corona_Cases.long,case_type=="death"),c("Country.Region","Province.State","City",all_of(metadata_columns))))
Corona_Cases.long<-select(Corona_Cases.long,-all_of(metadata_columns))

# Some counties are not summarized on the country level. collapse all but US
Corona_Cases.long<-rbind.fill(ddply(filter(Corona_Cases.long,!Country.Region=="US_state"),c("case_type","Country.Region","Date","Date.numeric"),summarise,Total_confirmed_cases=sum(Total_confirmed_cases)),filter(Corona_Cases.long,Country.Region=="US_state"))

# Put total case and deaths side-by-side (wide)
Corona_Cases<-spread(Corona_Cases.long,key = case_type,value = Total_confirmed_cases)

#Compute moratlity rate
Corona_Cases$mortality_rate<-Corona_Cases$death/Corona_Cases$total

#TMP
Corona_Cases<-plyr::rename(Corona_Cases,c("total"="Total_confirmed_cases","death"="Total_confirmed_deaths"))

##------------------------------------------
## log10 transform total # cases
##------------------------------------------
Corona_Cases$Total_confirmed_cases.log<-log(Corona_Cases$Total_confirmed_cases,10)
Corona_Cases$Total_confirmed_deaths.log<-log(Corona_Cases$Total_confirmed_deaths,10)
##------------------------------------------
       
##------------------------------------------
## Compute # of days since 100th for US data
##------------------------------------------

# Find day that 100th case was found for Country/Province. NOTE: Non US countries may have weird provinces. For example, Frane is summairzed at the country level but also had 3 providences. I've only ensured the U.S. case100 works... so the case100_date for U.S. is summarized both for the entire country (regardless of state) and on a per-state level. 
# TODO: consider city-level summary as well. This data may be sparse

Corona_Cases<-merge(Corona_Cases,ddply(filter(Corona_Cases,Total_confirmed_cases>100),c("Country.Region"),summarise,case100_date=min(Date.numeric)))
Corona_Cases$Days_since_100<-Corona_Cases$Date.numeric-Corona_Cases$case100_date

##------------------------------------------
## Add population and lat/long data (CURRENTLY US ONLY)
##------------------------------------------

kable(filter(metadata,(is.na(Country.Region) | is.na(Population) )) %>% select(c("Country.Region","Province.State","City")) %>% unique(),caption = "Regions for which either population or Country is NA")
Regions for which either population or Country is NA
Country.Region Province.State City
# Drop missing data 
metadata<-filter(metadata,!(is.na(Country.Region) | is.na(Population) ))
# Convert remaining pop to numeric
metadata$Population<-as.numeric(metadata$Population)
## Warning: NAs introduced by coercion
# Add metadata to cases
Corona_Cases<-merge(Corona_Cases,metadata,all.x = T)

##------------------------------------------
## Compute total and death cases relative to population 
##------------------------------------------

Corona_Cases$Total_confirmed_cases.per100<-100*Corona_Cases$Total_confirmed_cases/Corona_Cases$Population
Corona_Cases$Total_confirmed_deaths.per100<-100*Corona_Cases$Total_confirmed_deaths/Corona_Cases$Population


##------------------------------------------
## Filter df for US state-wide stats
##------------------------------------------

Corona_Cases.US_state<-filter(Corona_Cases,Country.Region=="US_state" & Total_confirmed_cases>0 ) 
kable(table(select(Corona_Cases.US_state,c("Province.State"))),caption = "Number of longitudinal datapoints (total/death) per state")
Number of longitudinal datapoints (total/death) per state
Var1 Freq
Alabama 4065
Alaska 747
Arizona 1031
Arkansas 4230
California 3933
Colorado 3706
Connecticut 610
Delaware 251
Diamond Princess 68
District of Columbia 69
Florida 4362
Georgia 9688
Grand Princess 69
Guam 69
Hawaii 357
Idaho 1925
Illinois 5502
Indiana 5641
Iowa 5085
Kansas 4173
Kentucky 6119
Louisiana 4085
Maine 1035
Maryland 1601
Massachusetts 1043
Michigan 4915
Minnesota 4588
Mississippi 5062
Missouri 5538
Montana 1730
Nebraska 3114
Nevada 765
New Hampshire 719
New Jersey 1553
New Mexico 1684
New York 3917
North Carolina 5936
North Dakota 1966
Northern Mariana Islands 54
Ohio 5301
Oklahoma 3976
Oregon 2042
Pennsylvania 4199
Puerto Rico 69
Rhode Island 406
South Carolina 2961
South Dakota 2509
Tennessee 5719
Texas 11862
Utah 954
Vermont 960
Virgin Islands 69
Virginia 7462
Washington 2683
West Virginia 2694
Wisconsin 3986
Wyoming 1231
Corona_Cases.US_state<-merge(Corona_Cases.US_state,ddply(filter(Corona_Cases.US_state,Total_confirmed_cases>100),c("Province.State"),summarise,case100_date_state=min(Date.numeric)))
Corona_Cases.US_state$Days_since_100_state<-Corona_Cases.US_state$Date.numeric-Corona_Cases.US_state$case100_date_state

ANALYSIS

Q1: What is the trend in cases, mortality across geopgraphical regions?

Plot # of cases vs time
* For each geographical set:
* comparative longitudinal case trend (absolute & log scale)
* comparative longitudinal mortality trend
* death vs total correlation

question dataset x y color facet pch dimentions
comparative_longitudinal_case_trend long time log_cases geography none (case type?) case_type [15, 50, 4] geography x (2 scale?) case type
comparative longitudinal case trend long time cases geography case_type ? [15, 50, 4] geography x (2+ scale) case type
comparative longitudinal mortality trend wide time mortality rate geography none none [15, 50, 4] geography
death vs total correlation wide cases deaths geography none none [15, 50, 4] geography
# total cases vs time
# death cases vs time
# mortality rate vs time
# death vs mortality


  # death vs mortality
  # total & death case vs time (same plot)

#<question> <x> <y> <colored> <facet> <dataset>
## trend in case/deaths over time, comapred across regions <time> <log cases> <geography*> <none> <.wide>
## trend in case/deaths over time, comapred across regions <time> <cases> <geography*> <case_type> <.long>
## trend in mortality rate over time, comapred across regions <time> <mortality rate> <geography*> <none>
## how are death/mortality related/correlated? <time> <log cases> <geography*> <none>
## how are death and case load correlated? <cases> <deaths>

# lm for each?? - > apply lm from each region starting from 100th case. m, b associated with each.
    # input: geographical regsion, logcase vs day (100th case)
    # output: m, b for each geographical region ID



#total/death on same plot-  diffeer by 2 logs, so when plotting log, use pch. when plotting absolute, need to use free scales
#when plotting death and case on same, melt. 

#CoronaCases - > filter sets (3)
  #world - choose countries with sufficent data

N<-ddply(filter(Corona_Cases,Total_confirmed_cases>100),c("Country.Region"),summarise,n=length(Country.Region))
ggplot(filter(N,n<100),aes(x=n))+
  geom_histogram()+
  default_theme+
  ggtitle("Distribution of number of days with at least 100 confirmed cases for each region")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

kable(arrange(N,-n),caption="Sorted number of days with at least 100 confirmed cases")
Sorted number of days with at least 100 confirmed cases
Country.Region n
US_state 35562
China 123
Diamond Princess 104
Korea, South 94
Japan 93
Italy 91
Iran 88
Singapore 85
France 84
Germany 84
Spain 83
US 82
Switzerland 80
United Kingdom 80
Belgium 79
Netherlands 79
Norway 79
Sweden 79
Austria 77
Malaysia 76
Australia 75
Bahrain 75
Denmark 75
Canada 74
Qatar 74
Iceland 73
Brazil 72
Czechia 72
Finland 72
Greece 72
Iraq 72
Israel 72
Portugal 72
Slovenia 72
Egypt 71
Estonia 71
India 71
Ireland 71
Kuwait 71
Philippines 71
Poland 71
Romania 71
Saudi Arabia 71
Indonesia 70
Lebanon 70
San Marino 70
Thailand 70
Chile 69
Pakistan 69
Luxembourg 68
Peru 68
Russia 68
Ecuador 67
Mexico 67
Slovakia 67
South Africa 67
United Arab Emirates 67
Armenia 66
Colombia 66
Croatia 66
Panama 66
Serbia 66
Taiwan* 66
Turkey 66
Argentina 65
Bulgaria 65
Latvia 65
Uruguay 65
Algeria 64
Costa Rica 64
Dominican Republic 64
Hungary 64
Andorra 63
Bosnia and Herzegovina 63
Jordan 63
Lithuania 63
Morocco 63
New Zealand 63
North Macedonia 63
Vietnam 63
Albania 62
Cyprus 62
Malta 62
Moldova 62
Brunei 61
Burkina Faso 61
Sri Lanka 61
Tunisia 61
Ukraine 60
Azerbaijan 59
Ghana 59
Kazakhstan 59
Oman 59
Senegal 59
Venezuela 59
Afghanistan 58
Cote d’Ivoire 58
Cuba 57
Mauritius 57
Uzbekistan 57
Cambodia 56
Cameroon 56
Honduras 56
Nigeria 56
West Bank and Gaza 56
Belarus 55
Georgia 55
Bolivia 54
Kosovo 54
Kyrgyzstan 54
Montenegro 54
Congo (Kinshasa) 53
Kenya 52
Niger 51
Guinea 50
Rwanda 50
Trinidad and Tobago 50
Paraguay 49
Bangladesh 48
Djibouti 46
El Salvador 45
Guatemala 44
Madagascar 43
Mali 42
Congo (Brazzaville) 39
Jamaica 39
Gabon 37
Somalia 37
Tanzania 37
Ethiopia 36
Burma 35
Sudan 34
Liberia 33
Maldives 31
Equatorial Guinea 30
Cabo Verde 28
Sierra Leone 26
Guinea-Bissau 25
Togo 25
Zambia 24
Eswatini 23
Chad 22
Tajikistan 21
Haiti 19
Sao Tome and Principe 19
Benin 17
Nepal 17
Uganda 17
Central African Republic 16
South Sudan 16
Guyana 14
Mozambique 13
Yemen 9
Mongolia 8
Mauritania 5
Nicaragua 5
# Pick top 15 countries with data
max_colors<-12
# find way to fix this- China has diff provences. Plot doesnt look right...
sufficient_data<-arrange(filter(N,!Country.Region %in% c("US_state", "Diamond Princess")),-n)[1:max_colors,]
kable(sufficient_data,caption = paste0("Top ",max_colors," countries with sufficient data"))
Top 12 countries with sufficient data
Country.Region n
China 123
Korea, South 94
Japan 93
Italy 91
Iran 88
Singapore 85
France 84
Germany 84
Spain 83
US 82
Switzerland 80
United Kingdom 80
Corona_Cases.world<-filter(Corona_Cases,Country.Region %in% c(sufficient_data$Country.Region))


  #us 
  #    - by state
Corona_Cases.US<-filter(Corona_Cases,Country.Region=="US" & Total_confirmed_cases>0)
# summarize 
#!City %in% c("Unassigned") 
  #    - specific cities
#mortality_rate!=Inf & mortality_rate<=1
Corona_Cases.UScity<-filter(Corona_Cases,Province.State %in% c("Pennsylvania","Maryland","New York","New Jersey") & City %in% c("Bucks","Baltimore City", "New York","Burlington","Cape May"))

measure_vars_long<-c("Total_confirmed_cases.log","Total_confirmed_cases","Total_confirmed_deaths","Total_confirmed_deaths.log")
melt_arg_list<-list(variable.name = "case_type",value.name = "cases",measure.vars = c("Total_confirmed_cases","Total_confirmed_deaths"))
melt_arg_list$data=NULL


melt_arg_list$data=select(Corona_Cases.world,-ends_with(match = "log"))
Corona_Cases.world.long<-do.call(melt,melt_arg_list)
melt_arg_list$data=select(Corona_Cases.UScity,-ends_with(match = "log"))
Corona_Cases.UScity.long<-do.call(melt,melt_arg_list)
melt_arg_list$data=select(Corona_Cases.US_state,-ends_with(match = "log"))
Corona_Cases.US_state.long<-do.call(melt,melt_arg_list)

Corona_Cases.world.long$cases.log<-log(Corona_Cases.world.long$cases,10)
Corona_Cases.US_state.long$cases.log<-log(Corona_Cases.US_state.long$cases,10)
Corona_Cases.UScity.long$cases.log<-log(Corona_Cases.UScity.long$cases,10)


# what is the current death and total case load for US? For world? For states?
#-absolute
#-log

# what is mortality rate (US, world)
#-absolute

#how is death and case correlated? (US, world)
#-absolute
#Corona_Cases.US<-filter(Corona_Cases,Country.Region=="US" & Total_confirmed_cases>0)
#Corona_Cases.US.case100<-filter(Corona_Cases.US, Days_since_100>=0)
# linear model parameters
#(model_fit<-lm(formula = Total_confirmed_cases.log~Days_since_100,data= Corona_Cases.US.case100 ))

#(slope<-model_fit$coefficients[2])
#(intercept<-model_fit$coefficients[1])

# Correlation coefficient
#cor(x = Corona_Cases.US.case100$Days_since_100,y = Corona_Cases.US.case100$Total_confirmed_cases.log)

##------------------------------------------
## Plot World Data
##------------------------------------------
# Timestamp for world
timestamp_plot.world<-paste("Most recent date for which data available:",max(Corona_Cases.world$Date))#timestamp(quiet = T,prefix = "Updated ",suffix = " (EST)")


# Base template for plots
baseplot.world<-ggplot(data=NULL,aes(x=Days_since_100,col=Country.Region))+
  default_theme+
  scale_color_brewer(type = "qualitative",palette = "Paired")+
  ggtitle(paste("Log10 cases over time,",timestamp_plot.world))+
  theme(legend.position = "bottom",plot.title = element_text(size=12))


##/////////////////////////
### Plot Longitudinal cases

(Corona_Cases.world.long.plot<-baseplot.world+
    geom_point(data=Corona_Cases.world.long,aes(y=cases))+
    geom_line(data=Corona_Cases.world.long,aes(y=cases))+
    facet_wrap(~case_type,scales = "free_y",ncol=1)+
    ggtitle(timestamp_plot.world)
    )

(Corona_Cases.world.loglong.plot<-baseplot.world+
    geom_point(data=Corona_Cases.world.long,aes(y=cases.log))+
    geom_line(data=Corona_Cases.world.long,aes(y=cases.log))+
    facet_wrap(~case_type,scales = "free_y",ncol=1)+
    ggtitle(timestamp_plot.world))

##/////////////////////////
### Plot Longitudinal mortality rate

(Corona_Cases.world.mortality.plot<-baseplot.world+
    geom_point(data=Corona_Cases.world,aes(y=mortality_rate))+
    geom_line(data=Corona_Cases.world,aes(y=mortality_rate))+
    ylim(c(0,0.3))+
    ggtitle(timestamp_plot.world))
## Warning: Removed 100 rows containing missing values (geom_point).
## Warning: Removed 100 row(s) containing missing values (geom_path).

##/////////////////////////
### Plot death vs total case correlation

(Corona_Cases.world.casecor.plot<-ggplot(Corona_Cases.world,aes(x=Total_confirmed_cases,y=Total_confirmed_deaths,col=Country.Region))+
  geom_point()+
  geom_line()+
  default_theme+
  scale_color_brewer(type = "qualitative",palette = "Paired")+
  ggtitle(paste("Log10 cases over time,",timestamp_plot.world))+
  theme(legend.position = "bottom",plot.title = element_text(size=12))+
    ggtitle(timestamp_plot.world))

### Write polots

write_plot(Corona_Cases.world.long.plot,wd = results_dir)
## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.long.plot.png"
write_plot(Corona_Cases.world.loglong.plot,wd = results_dir)
## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.loglong.plot.png"
write_plot(Corona_Cases.world.mortality.plot,wd = results_dir)
## Warning: Removed 100 rows containing missing values (geom_point).

## Warning: Removed 100 row(s) containing missing values (geom_path).
## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.mortality.plot.png"
write_plot(Corona_Cases.world.casecor.plot,wd = results_dir)
## [1] "/Users/stevensmith/Projects/coronavirus/results/Corona_Cases.world.casecor.plot.png"
##------------------------------------------
## Plot US State Data
##-----------------------------------------

baseplot.US<-ggplot(data=NULL,aes(x=Days_since_100_state,col=case_type))+
  default_theme+
  facet_wrap(~Province.State)+
  ggtitle(paste("Log10 cases over time,",timestamp_plot.world))

Corona_Cases.US_state.long.plot<-baseplot.US+geom_point(data=Corona_Cases.US_state.long,aes(y=cases.log))
##------------------------------------------
## Plot US City Data
##-----------------------------------------

Corona_Cases.US.plotdata<-filter(Corona_Cases.US_state,Province.State %in% c("Pennsylvania","Maryland","New York","New Jersey") &
                                   City %in% c("Bucks","Baltimore City", "New York","Burlington","Cape May") &
                                   Total_confirmed_cases>0) 
timestamp_plot<-paste("Most recent date for which data available:",max(Corona_Cases.US.plotdata$Date))#timestamp(quiet = T,prefix = "Updated ",suffix = " (EST)")

city_colors<-c("Bucks"='#beaed4',"Baltimore City"='#386cb0', "New York"='#7fc97f',"Burlington"='#fdc086',"Cape May"="#e78ac3")

##/////////////////////////
### Plot death vs total case correlation

(Corona_Cases.city.loglong.plot<-ggplot(melt(Corona_Cases.US.plotdata,measure.vars = c("Total_confirmed_cases.log","Total_confirmed_deaths.log"),variable.name = "case_type",value.name = "cases"),aes(x=Date,y=cases,col=City,pch=case_type))+
  geom_point(size=4)+
    geom_line()+
  default_theme+
  #facet_wrap(~case_type)+
    ggtitle(paste("Log10 total and death cases over time,",timestamp_plot))+
theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))+
    scale_color_manual(values = city_colors)+
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

(Corona_Cases.city.long.plot<-ggplot(filter(Corona_Cases.US.plotdata,Province.State !="New York"),aes(x=Date,y=Total_confirmed_cases,col=City))+
  geom_point(size=4)+
  geom_line()+
  default_theme+
  facet_grid(~Province.State,scales = "free_y")+
  ggtitle(paste("MD, PA, NJ total cases over time,",timestamp_plot))+
  theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))
+
  scale_color_manual(values = city_colors)+
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

(Corona_Cases.city.mortality.plot<-ggplot(Corona_Cases.US.plotdata,aes(x=Date,y=mortality_rate,col=City))+
  geom_point(size=3)+
  geom_line(size=2)+
  default_theme+
  ggtitle(paste("Mortality rate (deaths/total) over time,",timestamp_plot))+
  theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))+
  scale_color_manual(values = city_colors)+
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

(Corona_Cases.city.casecor.plot<-ggplot(filter(Corona_Cases.US.plotdata,Province.State !="New York"),aes(y=Total_confirmed_deaths,x=Total_confirmed_cases,col=City))+
  geom_point(size=3)+
  geom_line(size=2)+
  default_theme+
  ggtitle(paste("Correlation of death vs total cases,",timestamp_plot))+
  theme(legend.position = "bottom",plot.title = element_text(size=12))+
  scale_color_manual(values = city_colors))

(Corona_Cases.city.long.normalized.plot<-ggplot(filter(Corona_Cases.US.plotdata,Province.State !="New York"),aes(x=Date,y=Total_confirmed_cases.per100,col=City))+
  geom_point(size=4)+
  geom_line()+
  default_theme+
  facet_grid(~Province.State)+
  ggtitle(paste("MD, PA, NJ total cases over time per 100 people,",timestamp_plot))+
  theme(legend.position = "bottom",plot.title = element_text(size=12),axis.text.x = element_text(angle=45,hjust=1))+
  scale_color_manual(values = city_colors)  +
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

write_plot(Corona_Cases.city.long.plot,wd = results_dir_custom)
## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.long.plot.png"
write_plot(Corona_Cases.city.loglong.plot,wd = results_dir_custom)
## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.loglong.plot.png"
write_plot(Corona_Cases.city.mortality.plot,wd = results_dir_custom)
## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.mortality.plot.png"
write_plot(Corona_Cases.city.casecor.plot,wd = results_dir_custom)
## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.casecor.plot.png"
write_plot(Corona_Cases.city.long.normalized.plot,wd = results_dir_custom)
## [1] "/Users/stevensmith/Projects/coronavirus/results/custom/Corona_Cases.city.long.normalized.plot.png"

Q1b what is the model

Fit the cases to a linear model 1. Find time at which the case vs date becomes linear in each plot
2. Fit linear model for each city

# What is the predict # of cases for the next few days?
# How is the model performing historically?

Corona_Cases.US_state.summary<-ddply(Corona_Cases.US_state,
                                     c("Province.State","Date"),
                                     summarise,
                                     Total_confirmed_cases_perstate=sum(Total_confirmed_cases)) %>% 
    filter(Total_confirmed_cases_perstate>100)

# Compute the states with the most cases (for coloring and for linear model)
top_states_totals<-head(ddply(Corona_Cases.US_state.summary,c("Province.State"),summarise, Total_confirmed_cases_perstate.max=max(Total_confirmed_cases_perstate)) %>% arrange(-Total_confirmed_cases_perstate.max),n=max_colors)

kable(top_states_totals,caption = "Top 12 States, total count ")
top_states<-top_states_totals$Province.State

# Manually fix states so that Maryland is switched out for New York
top_states_modified<-c(top_states[top_states !="New York"],"Maryland")

# Plot with all states:
(Corona_Cases.US_state.summary.plot<-ggplot(Corona_Cases.US_state.summary,aes(x=Date,y=Total_confirmed_cases_perstate))+
  geom_point()+
  geom_point(data=filter(Corona_Cases.US_state.summary,Province.State %in% top_states),aes(col=Province.State))+
  scale_color_brewer(type = "qualitative",palette = "Paired")+
  default_theme+
  theme(axis.text.x = element_text(angle=45,hjust=1),legend.position = "bottom")+
  ggtitle("Total confirmed cases per state, top 12 colored")+
  scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

##------------------------------------------
## Fit linear model to time vs total cases
##-----------------------------------------

# First, find the date at which each state's cases vs time becomes lienar (2nd derivative is about 0)
li<-ddply(Corona_Cases.US_state.summary,c("Province.State"),find_linear_index)

# Compute linear model for each state starting at the point at which data becomes linear
for(i in 1:nrow(li)){
  Province.State.i<-li[i,"Province.State"]
  date.i<-li[i,"V1"]
  data.i<-filter(Corona_Cases.US_state.summary,Province.State==Province.State.i & as.numeric(Date) >= date.i)
  model_results<-lm(data.i,formula = Total_confirmed_cases_perstate~Date)
  slope<-model_results$coefficients[2]
  intercept<-model_results$coefficients[1]
  li[li$Province.State==Province.State.i,"m"]<-slope
  li[li$Province.State==Province.State.i,"b"]<-intercept
  }

# Compute top state case load with fitted model

(Corona_Cases.US_state.lm.plot<-ggplot(filter(Corona_Cases.US_state.summary,Province.State %in% top_states_modified ))+
    geom_abline(data=filter(li,Province.State %in% top_states_modified),
                aes(slope = m,intercept = b,col=Province.State),lty=2)+
    geom_point(aes(x=Date,y=Total_confirmed_cases_perstate,col=Province.State))+
    scale_color_brewer(type = "qualitative",palette = "Paired")+
    default_theme+
    theme(axis.text.x = element_text(angle=45,hjust=1),legend.position = "bottom")+
    ggtitle("Total confirmed cases per state, top 12 colored")+
    scale_x_date(date_breaks="1 week",date_minor_breaks="1 day"))

##------------------------------------------
## Predict the number of total cases over the next week
##-----------------------------------------

predicted_days<-c(0,1,2,3,7)+as.numeric(as.Date("2020-04-20"))

predicted_days_df<-data.frame(matrix(ncol=3))
names(predicted_days_df)<-c("Province.State","days","Total_confirmed_cases_perstate")

# USe model parameters to estiamte case loads
for(state.i in top_states_modified){
  predicted_days_df<-rbind(predicted_days_df,
                           data.frame(Province.State=state.i,
                                      prediction_model(m = li[li$Province.State==state.i,"m"],
                                                       b =li[li$Province.State==state.i,"b"] ,
                                                       days =predicted_days )))
  }

predicted_days_df$Date<-as.Date(predicted_days_df$days,origin="1970-01-01")

kable(predicted_days_df,caption = "Predicted total cases over the next week for selected states")

##------------------------------------------
## Write plots
##-----------------------------------------

write_plot(Corona_Cases.US_state.summary.plot,wd = results_dir)
write_plot(Corona_Cases.US_state.lm.plot,wd = results_dir)

##------------------------------------------
## Write tables
##-----------------------------------------

write.csv(predicted_days_df,file = paste0(results_dir,"predicted_total_cases_days.csv"),quote = F,row.names = F)

Q2: What is the predicted number of cases?

What is the prediction of COVID-19 based on model thus far? Additional questions:

WHy did it take to day 40 to start a log linear trend? How long will it be till x number of cases? When will the plateu happen? Are any effects noticed with social distancing? Delays

##------------------------------------------
## Prediction and Prediction Accuracy
##------------------------------------------


today_num<-max(Corona_Cases.US$Days_since_100)
predicted_days<-today_num+c(1,2,3,7)

#mods = dlply(mydf, .(x3), lm, formula = y ~ x1 + x2)
#today:
Corona_Cases.US[Corona_Cases.US$Days_since_100==(today_num-1),]
Corona_Cases.US[Corona_Cases.US$Days_since_100==today_num,]
Corona_Cases.US$type<-"Historical"


#prediction_values<-prediction_model(m=slope,b=intercept,days = predicted_days)$Total_confirmed_cases

histoical_model<-data.frame(date=today_num,m=slope,b=intercept)
tmp<-data.frame(state=rep(c("A","B"),each=3),x=c(1,2,3,4,5,6))
tmp$y<-c(tmp[1:3,"x"]+5,tmp[4:6,"x"]*5+1)
ddply(tmp,c("state"))
lm(data =tmp,formula = y~x )

train_lm<-function(input_data,subset_coulmn,formula_input){
case_models <- dlply(input_data, subset_coulmn, lm, formula = formula_input)
case_models.parameters <- ldply(case_models, coef)
case_models.parameters<-rename(case_models.parameters,c("b"="(Intercept)","m"=subset_coulmn))
return(case_models.parameters)
}

train_lm(tmp,"state")

 dlply(input_data, subset_coulmn, lm,m=)
 
# model for previous y days
#historical_model_predictions<-data.frame(day_x=NULL,Days_since_100=NULL,Total_confirmed_cases=NULL,Total_confirmed_cases.log=NULL)
# for(i in c(1,2,3,4,5,6,7,8,9,10)){
#   #i<-1
# day_x<-today_num-i # 1, 2, 3, 4
# day_x_nextweek<-day_x+c(1,2,3)
# model_fit_x<-lm(data = filter(Corona_Cases.US.case100,Days_since_100 < day_x),formula = Total_confirmed_cases.log~Days_since_100)
# prediction_day_x_nextweek<-prediction_model(m = model_fit_x$coefficients[2],b = model_fit_x$coefficients[1],days = day_x_nextweek)
# prediction_day_x_nextweek$type<-"Predicted"
# acutal_day_x_nextweek<-filter(Corona_Cases.US,Days_since_100 %in% day_x_nextweek) %>% select(c(Days_since_100,Total_confirmed_cases,Total_confirmed_cases.log))
# acutal_day_x_nextweek$type<-"Historical"
# historical_model_predictions.i<-data.frame(day_x=day_x,rbind(acutal_day_x_nextweek,prediction_day_x_nextweek))
# historical_model_predictions<-rbind(historical_model_predictions.i,historical_model_predictions)
# }

#historical_model_predictions.withHx<-rbind.fill(historical_model_predictions,data.frame(Corona_Cases.US,type="Historical"))
#historical_model_predictions.withHx$Total_confirmed_cases.log2<-log(historical_model_predictions.withHx$Total_confirmed_cases,2)

(historical_model_predictions.plot<-ggplot(historical_model_predictions.withHx,aes(x=Days_since_100,y=Total_confirmed_cases.log,col=type))+
    geom_point(size=3)+
    default_theme+
    theme(legend.position = "bottom")+ 
      #geom_abline(slope = slope,intercept =intercept,lty=2)+
    #facet_wrap(~case_type,ncol=1)+
    scale_color_manual(values = c("Historical"="#377eb8","Predicted"="#e41a1c")))
write_plot(historical_model_predictions.plot,wd=results_dir)

Q3: What is the effect on social distancing, descreased mobility on case load?

Load data from Google which compoutes % change in user mobility relative to baseline for * Recreation
* Workplace
* Residence
* Park
* Grocery

Data from https://www.google.com/covid19/mobility/

# See pre-processing section for script on gathering mobility data

# UNDER DEVELOPMENT

mobility<-read.csv("/Users/stevensmith/Projects/MIT_COVID19/mobility.csv",header = T,stringsAsFactors = F)
#mobility$Retail_Recreation<-as.numeric(sub(mobility$Retail_Recreation,pattern = "%",replacement = ""))
#mobility$Workplace<-as.numeric(sub(mobility$Workplace,pattern = "%",replacement = ""))
#mobility$Residential<-as.numeric(sub(mobility$Residential,pattern = "%",replacement = ""))

##------------------------------------------
## Show relationship between mobility and caseload
##------------------------------------------
mobility$County<-gsub(mobility$County,pattern = " County",replacement = "")
Corona_Cases.US_state.mobility<-merge(Corona_Cases.US_state,plyr::rename(mobility,c("State"="Province.State","County"="City")))

#Corona_Cases.US_state.tmp<-merge(metadata,Corona_Cases.US_state.tmp)
# Needs to happen upsteam, see todos
#Corona_Cases.US_state.tmp$Total_confirmed_cases.perperson<-Corona_Cases.US_state.tmp$Total_confirmed_cases/as.numeric(Corona_Cases.US_state.tmp$Population)
mobility_measures<-c("Retail_Recreation","Grocery_Pharmacy","Parks","Transit","Workplace","Residential")

plot_data<-filter(Corona_Cases.US_state.mobility, Date.numeric==max(Corona_Cases.US_state$Date.numeric) ) %>% melt(measure.vars=mobility_measures) 
plot_data$value<-as.numeric(gsub(plot_data$value,pattern = "%",replacement = ""))
plot_data<-filter(plot_data,!is.na(value))

(mobility.plot<-ggplot(filter(plot_data,Province.State %in% c("Pennsylvania","Maryland","New Jersey","California","Delaware","Connecticut")),aes(y=Total_confirmed_cases.per100,x=value))+geom_point()+
  facet_grid(Province.State~variable,scales = "free")+
  xlab("Mobility change from baseline (%)")+
  ylab(paste0("Confirmed cases per 100 people(Today)"))+
  default_theme+
  ggtitle("Mobility change vs cases"))

(mobility.global.plot<-ggplot(plot_data,aes(y=Total_confirmed_cases.per100,x=value))+geom_point()+
  facet_wrap(~variable,scales = "free")+
  xlab("Mobility change from baseline (%)")+
  ylab(paste0("Confirmed cases (Today) per 100 people"))+
  default_theme+
  ggtitle("Mobility change vs cases"))

plot_data.permobility_summary<-ddply(plot_data,c("Province.State","variable"),summarise,cor=cor(y =Total_confirmed_cases.per100,x=value),median_change=median(x=value)) %>% arrange(-abs(cor))

kable(plot_data.permobility_summary,caption = "Ranked per-state mobility correlation with total confirmed cases")
Ranked per-state mobility correlation with total confirmed cases
Province.State variable cor median_change
Alaska Transit -1.0000000 -63.0
Delaware Retail_Recreation 1.0000000 -39.5
Delaware Grocery_Pharmacy 1.0000000 -17.5
Delaware Parks -1.0000000 20.5
Delaware Transit 1.0000000 -37.0
Delaware Workplace 1.0000000 -37.0
Delaware Residential -1.0000000 14.0
Hawaii Retail_Recreation 0.9931972 -56.0
Hawaii Grocery_Pharmacy 0.9695437 -34.0
New Hampshire Parks 0.9569819 -20.0
Connecticut Grocery_Pharmacy -0.9034799 -6.0
Maine Transit -0.9024399 -50.0
Alaska Residential 0.8898278 13.0
South Dakota Parks 0.8710801 -26.0
Utah Residential -0.8524968 12.0
Vermont Parks 0.8519884 -35.5
Alaska Grocery_Pharmacy -0.8078182 -7.0
Hawaii Residential -0.7854909 19.0
Utah Transit -0.7794922 -18.0
Massachusetts Workplace -0.7711813 -39.0
Connecticut Transit -0.7510764 -50.0
Rhode Island Workplace -0.7500031 -39.5
Alaska Workplace -0.7296146 -34.0
Wyoming Transit -0.7287206 -17.0
Utah Parks -0.6985837 17.0
Wyoming Parks -0.6881358 -4.0
Hawaii Parks 0.6813458 -72.0
Utah Workplace -0.6591668 -37.0
Vermont Grocery_Pharmacy -0.6585379 -25.0
New York Workplace -0.6458510 -34.5
Maine Workplace -0.6449558 -30.0
Rhode Island Retail_Recreation -0.6292286 -45.0
Montana Workplace -0.6239388 -40.5
Arizona Grocery_Pharmacy -0.6233132 -15.0
Hawaii Transit 0.6188732 -89.0
Rhode Island Residential -0.6159667 18.5
New Jersey Workplace -0.6117849 -44.0
Nebraska Workplace 0.6068569 -32.5
New Jersey Parks -0.5966373 -6.0
New York Retail_Recreation -0.5871506 -46.0
North Dakota Retail_Recreation -0.5596690 -42.0
Hawaii Workplace 0.5396454 -46.0
Connecticut Residential 0.5358815 14.0
Massachusetts Retail_Recreation -0.5293009 -44.0
New York Parks 0.5268744 20.0
Connecticut Retail_Recreation -0.5165179 -45.0
New Jersey Retail_Recreation -0.5105718 -62.5
Maine Parks 0.5013907 -31.0
Connecticut Workplace -0.4988833 -39.0
Arizona Retail_Recreation -0.4981312 -42.5
North Dakota Parks 0.4973765 -34.0
Montana Parks -0.4913929 -58.0
New Jersey Grocery_Pharmacy -0.4858699 2.5
Nebraska Residential -0.4858344 14.0
New Mexico Grocery_Pharmacy -0.4793089 -11.0
Wyoming Workplace -0.4786493 -31.0
Iowa Parks -0.4777179 28.5
Rhode Island Parks 0.4731974 52.0
Montana Residential 0.4701424 14.0
New Mexico Parks 0.4548477 -31.5
Illinois Transit -0.4517787 -31.0
New Mexico Residential 0.4491212 13.5
Vermont Residential 0.4379791 11.5
Idaho Workplace -0.4374387 -29.0
California Transit -0.4371164 -42.0
Kentucky Parks -0.4360494 28.5
California Residential 0.4329382 14.0
Pennsylvania Workplace -0.4324064 -36.0
Wisconsin Transit -0.4318169 -23.5
Massachusetts Grocery_Pharmacy -0.4313347 -7.0
New Jersey Transit -0.4278509 -50.5
Kansas Parks 0.4187728 72.0
South Carolina Workplace 0.4180317 -30.0
Montana Retail_Recreation -0.4145442 -51.0
New Hampshire Residential -0.4124232 14.0
Arkansas Parks -0.4086800 -12.0
Idaho Transit -0.4084972 -30.0
Idaho Grocery_Pharmacy -0.4036292 -4.5
Alabama Workplace -0.3999512 -29.0
Maryland Workplace -0.3981259 -35.0
Arizona Residential 0.3943607 13.0
Maryland Grocery_Pharmacy -0.3897732 -10.0
Montana Transit -0.3892824 -41.0
Alabama Grocery_Pharmacy -0.3814724 -2.0
New York Transit -0.3745865 -48.0
Arizona Transit 0.3632915 -38.0
Florida Residential 0.3592640 14.0
New Mexico Retail_Recreation -0.3574663 -42.5
Pennsylvania Retail_Recreation -0.3562258 -45.0
Wyoming Grocery_Pharmacy -0.3534260 -10.0
California Parks -0.3461554 -38.5
Alabama Transit -0.3436527 -36.5
Nevada Transit -0.3424973 -20.0
Michigan Parks 0.3370423 30.0
Nebraska Grocery_Pharmacy 0.3317520 -0.5
Pennsylvania Parks 0.3292215 13.0
Montana Grocery_Pharmacy -0.3279939 -16.0
Alaska Retail_Recreation 0.3260634 -39.0
West Virginia Parks 0.3239758 -33.0
Minnesota Transit -0.3157689 -28.5
North Carolina Grocery_Pharmacy 0.3100644 0.0
Idaho Retail_Recreation -0.3048434 -40.5
West Virginia Grocery_Pharmacy -0.3040290 -6.0
Maine Retail_Recreation -0.3010249 -42.0
Vermont Retail_Recreation 0.2952799 -57.0
North Dakota Workplace 0.2903163 -40.0
Colorado Residential 0.2883815 14.0
Texas Workplace 0.2807795 -32.0
California Retail_Recreation -0.2801884 -44.0
Utah Retail_Recreation -0.2800449 -40.0
Mississippi Residential 0.2786968 13.0
Texas Residential -0.2786580 15.0
Arkansas Retail_Recreation -0.2768985 -30.0
Maryland Retail_Recreation -0.2758363 -39.0
Rhode Island Transit -0.2752347 -56.0
Kansas Workplace 0.2739561 -32.5
California Grocery_Pharmacy -0.2723915 -11.5
Maryland Residential 0.2682304 15.0
Oregon Grocery_Pharmacy 0.2651842 -7.0
Florida Parks -0.2650683 -43.0
North Carolina Workplace 0.2648515 -31.0
Nevada Residential 0.2634265 17.0
California Workplace -0.2620151 -36.0
Nevada Retail_Recreation -0.2612579 -43.0
Virginia Transit -0.2600154 -33.0
Rhode Island Grocery_Pharmacy 0.2565848 -7.5
Illinois Workplace -0.2525444 -31.0
Texas Parks 0.2522483 -42.0
Tennessee Workplace -0.2513115 -31.0
Vermont Workplace -0.2476498 -43.0
Tennessee Residential 0.2469408 11.5
Wisconsin Parks 0.2459968 51.5
Georgia Grocery_Pharmacy -0.2451624 -10.0
South Carolina Parks -0.2442901 -23.0
Illinois Parks 0.2393778 26.5
Pennsylvania Grocery_Pharmacy -0.2367309 -6.0
New York Grocery_Pharmacy -0.2318705 8.0
Arkansas Residential 0.2316479 12.0
Washington Workplace -0.2199559 -38.0
Michigan Workplace -0.2171009 -40.0
Missouri Residential -0.2140702 13.0
North Carolina Transit 0.2137332 -32.0
North Carolina Residential 0.2119263 13.0
New Jersey Residential 0.2104578 18.0
Missouri Workplace 0.2011787 -28.5
Kansas Grocery_Pharmacy -0.1995442 -14.0
Iowa Transit 0.1974885 -24.0
Oregon Residential 0.1928806 10.5
South Dakota Transit -0.1921556 -40.0
Georgia Workplace -0.1891422 -33.5
Illinois Residential 0.1870788 14.0
Mississippi Grocery_Pharmacy -0.1840467 -8.0
North Dakota Grocery_Pharmacy -0.1836997 -8.0
Idaho Residential -0.1806540 11.0
Colorado Parks -0.1791649 2.0
Wyoming Retail_Recreation -0.1771708 -39.0
Alabama Parks 0.1718818 -1.0
Texas Transit 0.1693244 -41.5
Georgia Retail_Recreation -0.1692230 -41.0
New Mexico Transit 0.1679544 -38.5
Virginia Grocery_Pharmacy -0.1645141 -8.0
Ohio Transit 0.1630571 -28.0
Wisconsin Residential -0.1629294 14.0
Georgia Residential -0.1618675 13.0
Virginia Residential 0.1604840 14.0
Virginia Parks 0.1596643 6.0
Oklahoma Residential 0.1567643 15.0
South Carolina Residential -0.1538989 12.0
South Dakota Retail_Recreation -0.1506499 -38.5
Washington Transit -0.1504574 -33.5
Connecticut Parks 0.1474430 43.0
Florida Retail_Recreation 0.1465621 -43.0
Washington Residential 0.1462162 13.0
Massachusetts Transit -0.1460918 -45.0
North Dakota Transit 0.1433902 -48.0
Michigan Retail_Recreation -0.1432653 -53.0
New Hampshire Retail_Recreation -0.1429439 -41.0
Indiana Retail_Recreation 0.1394149 -38.0
Indiana Residential 0.1387390 12.0
Massachusetts Parks 0.1379458 39.0
Maine Residential -0.1376179 11.0
Mississippi Transit -0.1373815 -38.5
North Carolina Parks -0.1315740 7.0
Florida Workplace -0.1314709 -33.0
Oregon Retail_Recreation 0.1305017 -41.0
Pennsylvania Transit -0.1302178 -41.5
Alabama Retail_Recreation 0.1295806 -39.0
South Dakota Residential 0.1255465 15.0
Oklahoma Parks -0.1252441 -18.5
Minnesota Parks 0.1245200 -9.0
Massachusetts Residential 0.1227471 15.0
Ohio Parks -0.1206226 67.5
New Hampshire Grocery_Pharmacy -0.1185951 -6.0
Wisconsin Workplace -0.1183345 -31.0
Kansas Transit -0.1168930 -26.5
Maine Grocery_Pharmacy -0.1167170 -13.0
Washington Grocery_Pharmacy 0.1163598 -7.0
Minnesota Workplace -0.1149393 -33.0
Oregon Parks 0.1143547 16.5
Idaho Parks 0.1110882 -22.0
Maryland Transit -0.1101439 -39.0
Mississippi Retail_Recreation -0.1099643 -40.0
Kentucky Grocery_Pharmacy 0.1098535 4.0
Arkansas Transit 0.1090221 -27.0
West Virginia Residential -0.1087555 11.0
Arkansas Workplace -0.1080883 -26.0
Ohio Residential 0.1077552 14.0
Indiana Parks -0.1023437 29.0
New Hampshire Transit -0.1019342 -57.0
Arizona Workplace -0.1006545 -35.0
Nebraska Retail_Recreation 0.0992677 -36.0
Minnesota Retail_Recreation 0.0973131 -40.0
Mississippi Workplace -0.0956580 -33.0
Wisconsin Grocery_Pharmacy 0.0954311 -1.0
Georgia Parks 0.0951224 -6.0
Michigan Grocery_Pharmacy -0.0931909 -11.0
South Dakota Workplace 0.0927870 -35.0
Texas Grocery_Pharmacy 0.0923568 -14.0
Wyoming Residential 0.0920486 12.5
Pennsylvania Residential 0.0915801 15.0
New York Residential 0.0913396 17.5
Nebraska Parks 0.0911200 55.5
Oklahoma Grocery_Pharmacy -0.0895337 -1.0
Missouri Transit -0.0881870 -24.5
Indiana Workplace 0.0871502 -34.0
Virginia Workplace -0.0869857 -31.5
Washington Parks 0.0859680 -3.5
South Dakota Grocery_Pharmacy 0.0843545 -9.0
Kentucky Transit 0.0794671 -31.0
Michigan Residential 0.0779151 15.0
Indiana Grocery_Pharmacy -0.0766021 -5.5
South Carolina Transit 0.0760657 -45.0
Kentucky Retail_Recreation 0.0747466 -29.0
Virginia Retail_Recreation -0.0742280 -35.0
Colorado Transit 0.0719691 -36.0
Ohio Grocery_Pharmacy 0.0666812 0.0
North Carolina Retail_Recreation 0.0662733 -34.0
Nevada Workplace 0.0655054 -40.0
Michigan Transit 0.0651504 -46.0
Tennessee Parks -0.0645235 10.5
Ohio Retail_Recreation 0.0622588 -36.0
West Virginia Workplace 0.0612386 -33.0
Minnesota Grocery_Pharmacy 0.0599051 -6.0
Oregon Transit 0.0597446 -27.5
Nevada Parks 0.0594394 -12.5
Oregon Workplace -0.0589402 -31.0
Nebraska Transit -0.0576652 -9.0
Iowa Retail_Recreation -0.0565372 -38.0
West Virginia Retail_Recreation -0.0526565 -38.5
Florida Transit -0.0516869 -49.0
Oklahoma Workplace 0.0512421 -31.0
South Carolina Retail_Recreation -0.0496295 -35.0
South Carolina Grocery_Pharmacy 0.0488516 1.0
Missouri Parks 0.0484492 0.0
Colorado Grocery_Pharmacy -0.0482725 -17.0
Colorado Retail_Recreation -0.0480970 -44.0
Washington Retail_Recreation -0.0480542 -42.0
Illinois Retail_Recreation 0.0475648 -40.0
Kentucky Residential 0.0473650 12.0
Texas Retail_Recreation 0.0470931 -40.0
New Hampshire Workplace 0.0452092 -37.0
Kentucky Workplace -0.0449154 -36.0
Missouri Grocery_Pharmacy 0.0441849 2.0
Missouri Retail_Recreation -0.0437999 -36.0
Iowa Workplace -0.0423004 -30.0
North Dakota Residential -0.0403637 17.0
Illinois Grocery_Pharmacy -0.0379889 2.0
Tennessee Transit -0.0376040 -32.0
Florida Grocery_Pharmacy 0.0335606 -14.0
Arizona Parks -0.0322416 -44.5
Indiana Transit 0.0291950 -29.0
Minnesota Residential -0.0288864 17.0
Vermont Transit 0.0270377 -63.0
Tennessee Grocery_Pharmacy 0.0262641 6.0
Ohio Workplace -0.0258992 -35.0
Oklahoma Retail_Recreation 0.0230294 -31.0
Nevada Grocery_Pharmacy 0.0189860 -12.5
West Virginia Transit -0.0170546 -45.0
Kansas Residential -0.0166925 13.0
Wisconsin Retail_Recreation 0.0156966 -44.0
Mississippi Parks -0.0145460 -25.0
Kansas Retail_Recreation -0.0131605 -37.0
Utah Grocery_Pharmacy 0.0130331 -4.0
Georgia Transit -0.0124298 -35.0
Maryland Parks -0.0105794 27.0
Iowa Residential -0.0081391 13.0
Colorado Workplace -0.0079635 -39.0
New Mexico Workplace 0.0071133 -34.0
Oklahoma Transit 0.0045994 -26.0
Tennessee Retail_Recreation -0.0036535 -30.0
Iowa Grocery_Pharmacy 0.0025521 4.0
Arkansas Grocery_Pharmacy 0.0024220 3.0
Alabama Residential -0.0009927 11.0
Alaska Parks NA 29.0
District of Columbia Retail_Recreation NA -69.0
District of Columbia Grocery_Pharmacy NA -28.0
District of Columbia Parks NA -65.0
District of Columbia Transit NA -69.0
District of Columbia Workplace NA -48.0
District of Columbia Residential NA 17.0
# sanity check
ggplot(filter(plot_data,Province.State %in% c("Pennsylvania","Maryland","New Jersey","California","Delaware","Connecticut")),aes(x=Total_confirmed_cases.per100,fill=variable))+geom_histogram()+
  facet_grid(~Province.State)+
    default_theme+
  theme(legend.position = "bottom")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

write_plot(mobility.plot,wd = results_dir)
## [1] "/Users/stevensmith/Projects/coronavirus/results/mobility.plot.png"
write_plot(mobility.global.plot,wd = results_dir)
## [1] "/Users/stevensmith/Projects/coronavirus/results/mobility.global.plot.png"
(plot_data.permobility_summary.plot<-ggplot(plot_data.permobility_summary,aes(x=variable,y=median_change))+
  geom_jitter(size=2,width=.2)+
  #geom_jitter(data=plot_data.permobility_summary %>% arrange(-abs(median_change)) %>% head(n=15),aes(col=Province.State),size=2,width=.2)+
  default_theme+
  ggtitle("Per-Sate Median Change in Mobility")+
  xlab("Mobility Meaure")+
  ylab("Median Change from Baseline"))

write_plot(plot_data.permobility_summary.plot,wd = results_dir)
## [1] "/Users/stevensmith/Projects/coronavirus/results/plot_data.permobility_summary.plot.png"

DELIVERABLE MANIFEST

The following link to commited documents pushed to github. These are provided as a convienence, but note this is a manual process. The generation of reports, plots and tables is not coupled to the execution of this markdown. ## Report This report, html & pdf

Plots

github_root<-"https://github.com/sbs87/coronavirus/blob/master/"

plot_handle<-c("Corona_Cases.world.long.plot",
               "Corona_Cases.world.loglong.plot",
               "Corona_Cases.world.mortality.plot",
               "Corona_Cases.world.casecor.plot",
               "Corona_Cases.city.long.plot",
               "Corona_Cases.city.loglong.plot",
               "Corona_Cases.city.mortality.plot",
               "Corona_Cases.city.casecor.plot",
               "Corona_Cases.city.long.normalized.plot",
               "Corona_Cases.US_state.lm.plot",
               "Corona_Cases.US_state.summary.plot")


deliverable_manifest<-data.frame(
  name=c("World total & death cases, longitudinal",
         "World log total & death cases, longitudinal",
         "World mortality",
         "World total & death cases, correlation",
         "City total & death cases, longitudinal",
         "City log total & death cases, longitudinal",
         "City mortality",
         "City total & death cases, correlation",
         "City population normalized total & death cases, longitudinal",
         "State total cases (select) with linear model, longitudinal",
         "State total cases, longitudinal"),
  plot_handle=plot_handle,
  link=paste0(github_root,"results/",plot_handle,".png")
)


(tmp<-data.frame(row_out=apply(deliverable_manifest,MARGIN = 1,FUN = function(x) paste(x[1],x[2],x[3],sep=" | "))))
##                                                                                                                                                                                                        row_out
## 1                                           World total & death cases, longitudinal | Corona_Cases.world.long.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.long.plot.png
## 2                                 World log total & death cases, longitudinal | Corona_Cases.world.loglong.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.loglong.plot.png
## 3                                                         World mortality | Corona_Cases.world.mortality.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.mortality.plot.png
## 4                                      World total & death cases, correlation | Corona_Cases.world.casecor.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.casecor.plot.png
## 5                                              City total & death cases, longitudinal | Corona_Cases.city.long.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.plot.png
## 6                                    City log total & death cases, longitudinal | Corona_Cases.city.loglong.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.loglong.plot.png
## 7                                                            City mortality | Corona_Cases.city.mortality.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.mortality.plot.png
## 8                                         City total & death cases, correlation | Corona_Cases.city.casecor.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.casecor.plot.png
## 9  City population normalized total & death cases, longitudinal | Corona_Cases.city.long.normalized.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.normalized.plot.png
## 10                     State total cases (select) with linear model, longitudinal | Corona_Cases.US_state.lm.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.lm.plot.png
## 11                                      State total cases, longitudinal | Corona_Cases.US_state.summary.plot | https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.summary.plot.png
row_out<-apply(tmp, 2, paste, collapse="\t\n")
name handle link
World total & death cases, longitudinal Corona_Cases.world.long.plot https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.long.plot.png
World log total & death cases, longitudinal Corona_Cases.world.loglong.plot https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.loglong.plot.png
World mortality Corona_Cases.world.mortality.plot https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.mortality.plot.png
World total & death cases, correlation Corona_Cases.world.casecor.plot https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.world.casecor.plot.png
City total & death cases, longitudinal Corona_Cases.city.long.plot https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.plot.png
City log total & death cases, longitudinal Corona_Cases.city.loglong.plot https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.loglong.plot.png
City mortality Corona_Cases.city.mortality.plot https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.mortality.plot.png
City total & death cases, correlation Corona_Cases.city.casecor.plot https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.casecor.plot.png
City population normalized total & death cases, longitudinal Corona_Cases.city.long.normalized.plot https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.city.long.normalized.plot.png
State total cases (select) with linear model, longitudinal Corona_Cases.US_state.lm.plot https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.lm.plot.png
State total cases, longitudinal Corona_Cases.US_state.summary.plot https://github.com/sbs87/coronavirus/blob/master/results/Corona_Cases.US_state.summary.plot.png

Tables

CONCLUSION

Overall, the trends of COVID-19 cases is no longer in log-linear phase for world or U.S. (but some regions like MD are still in the log-linear phase). Mortality rate (deaths/confirmed RNA-based cases) is >1%, with a range depending on region. Mobility is not a strong indicator of caseload (U.S. data).

See table below for detailed breakdown.

Question Answer
What is the effect on social distancing, descreased mobility on case load?
There is not a strong apparent effect on decreased mobility (work, grocery, retail) or increased mobility (at residence, parks) on number of confirmed cases, either as a country (U.S.) or state level. California appears to have one of the best correlations, but this is a mixed bag
What is the trend in cases, mortality across geopgraphical regions?
The confirmed total casees and mortality is overall log-linear for most countries, with a trailing off beginning for most (inlcuding U.S.). On the state level, NY, NJ, PA starting to trail off; MD is still in log-linear phase. Mortality and case load are highly correlated for NY, NJ, PA, MD. The mortality rate flucutates for a given region, but is about 3% overall.

END

End: ##—— Sun May 24 10:47:15 2020 ——##

Cheatsheet: http://rmarkdown.rstudio.com>

Sandbox

# Geographical heatmap!
install.packages("maps")
library(maps)
library
mi_counties <- map_data("county", "pennsylvania") %>% 
  select(lon = long, lat, group, id = subregion)
head(mi_counties)

ggplot(mi_counties, aes(lon, lat)) + 
  geom_point(size = .25, show.legend = FALSE) +
  coord_quickmap()
mi_counties$cases<-1:2226
name_overlaps(metadata,Corona_Cases.US_state)

tmp<-merge(Corona_Cases.US_state,metadata)
ggplot(filter(tmp,Province.State=="Pennsylvania"), aes(Long, Lat, group = as.factor(City))) +
  geom_polygon(aes(fill = Total_confirmed_cases), colour = "grey50") + 
  coord_quickmap()


ggplot(Corona_Cases.US_state, aes(Long, Lat))+
  geom_polygon(aes(fill = Total_confirmed_cases ), color = "white")+
  scale_fill_viridis_c(option = "C")
dev.off()


require(maps)
require(viridis)

world_map <- map_data("world")
ggplot(world_map, aes(x = long, y = lat, group = group)) +
  geom_polygon(fill="lightgray", colour = "white")

head(world_map)
head(Corona_Cases.US_state)
unique(select(world_map,c("region","group"))) %>% filter()

some.eu.countries <- c(
  "US"
)
# Retrievethe map data
some.eu.maps <- map_data("world", region = some.eu.countries)

# Compute the centroid as the mean longitude and lattitude
# Used as label coordinate for country's names
region.lab.data <- some.eu.maps %>%
  group_by(region) %>%
  summarise(long = mean(long), lat = mean(lat))

unique(filter(some.eu.maps,subregion %in% Corona_Cases.US_state$Province.State) %>% select(subregion))
unique(Corona_Cases.US_state$Total_confirmed_cases.log)
ggplot(filter(Corona_Cases.US_state,Date=="2020-04-17") aes(x = Long, y = Lat)) +
  geom_polygon(aes( fill = Total_confirmed_cases.log))+
  #geom_text(aes(label = region), data = region.lab.data,  size = 3, hjust = 0.5)+
  #scale_fill_viridis_d()+
  #theme_void()+
  theme(legend.position = "none")
library("sf")
library("rnaturalearth")
library("rnaturalearthdata")

world <- ne_countries(scale = "medium", returnclass = "sf")
class(world)
ggplot(data = world) +
    geom_sf()

counties <- st_as_sf(map("county", plot = FALSE, fill = TRUE))
counties <- subset(counties, grepl("florida", counties$ID))
counties$area <- as.numeric(st_area(counties))
#install.packages("lwgeom")
class(counties)
head(counties)
ggplot(data = world) +
    geom_sf(data=Corona_Cases.US_state) +
    #geom_sf(data = counties, aes(fill = area)) +
  geom_sf(data = counties, aes(fill = area)) +
   # scale_fill_viridis_c(trans = "sqrt", alpha = .4) +
    coord_sf(xlim = c(-88, -78), ylim = c(24.5, 33), expand = FALSE)


head(counties)
tmp<-unique(select(filter(Corona_Cases.US_state,Date=="2020-04-17"),c(Lat,Long,Total_confirmed_cases.per100)))
st_as_sf(map("county", plot = FALSE, fill = TRUE))

join::inner_join.sf(Corona_Cases.US_state, counties)

library(sf)
library(sp)

nc <- st_read(system.file("shape/nc.shp", package="sf"))
class(nc)


spdf <- SpatialPointsDataFrame(coords = select(Corona_Cases.US_state,c("Lat","Long")), data = Corona_Cases.US_state,
                               proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))

head(spdf)
class(spdf)
st_cast(spdf)

filter(Corona_Cases.US_state.summary,Date=="2020-04-20" & Province.State %in% top_states_modified)
id

https://stevenbsmith.net